JLA-ONA-2023

knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#install.packages("ona", repos = c("https://cran.qe-libs.org", "https://cran.rstudio.org")) # only run once
#install.packages("tma", repos = c("https://cran.qe-libs.org", "https://cran.rstudio.org")) # only run once

most recent data Ayano shared

data <- read.csv("~/Rprojects/ayano/PES_010v6(selected).csv")

units <- c("Speaker", "Role", "TimePhase")

codes_1 <- c("Functions", "Aesthetics", "User", "Vision", "Prototyping") # design actions codes
codes_2 <- c("CA", "ALoK", "CSU", "GCA", "Projective","Regulative","Relational") # shared epistemic agency codes

hoo_rules <- conversation_rules(
  (TimePhase %in% UNIT$TimePhase & Role %in% UNIT$Role)
)

notes: accum_1 and set_1 are for design action codes; accum_2 and set_2 are for shared epistemic agency codes

accum_1 <-
  contexts(data, 
           units_by = units, 
           hoo_rules = hoo_rules) %>%
  accumulate_contexts(codes = codes_1, 
                      decay.function = decay(simple_window, window_size = 4),
                      return.ena.set = FALSE, norm.by = NULL)

Note: Not doing means rotation so I can compare everything in the same space

set_1 <-
  # model(accum_1,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_1$meta.data$Role=="Engineer",
  #              servicedesigners=accum_1$meta.data$Role=="ServiceDesigner"))
  
  model(accum_1)

set_1 <-
  # model(accum_1,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_1$meta.data$Role=="Engineer",
  #              servicedesigners=accum_1$meta.data$Role=="ServiceDesigner"))
  
  model(accum_1)

set_2 <-
  # model(accum_2,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_2$meta.data$Role=="Engineer",
  #              servicedesigners=accum_2$meta.data$Role=="ServiceDesigner"))
  
  model(accum_2)

global visual parameter (let’s make sure all the plots use the same level of multiplier)

set_2 <-
  # model(accum_2,
  #       rotate.using = "mean",
  #       rotation.params = 
  #         list(engineers=accum_2$meta.data$Role=="Engineer",
  #              servicedesigners=accum_2$meta.data$Role=="ServiceDesigner"))
  
  model(accum_2)

Model 1: design action codes

Plot 1a: Overall Group Comparisons

traces = c(2:10)
plot(set_1, title = "Groups") |>
  units(
    points= set_1$points$Role$ProductDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("red"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_1$points$Role$Engineer, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("blue"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_1$points$Role$ServiceDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("green"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = "Product Designer", traces = traces[1]) |>
      plotly::style(name = "Engineer", traces = traces[2]) |>
      plotly::style(name = "Service Designer", traces = traces[3])

NA

Plot 1b: Group subtractions

prod_pts = set_1$points$Role$ProductDesigner 
eng_pts = set_1$points$Role$Engineer
serv_pts = set_1$points$Role$ServiceDesigner

prod_mean_net = set_1$line.weights %>% dplyr::filter(Role == "ProductDesigner") %>% colMeans()
eng_mean_net = set_1$line.weights %>% dplyr::filter(Role == "Engineer") %>% colMeans()
serv_mean_net = set_1$line.weights %>% dplyr::filter(Role == "ServiceDesigner") %>% colMeans()

plot(set_1, title = "Product Designer vs Engineer") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))

plot(set_1, title = "Product Designer vs ServiceDesigner") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points =serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - serv_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","green")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","green"))

plot(set_1, title = "Service Designer vs Engineer") |>
  units(
    points = serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = serv_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("green","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("green","blue"))

Plot 2: All Means by phase

#get list of time phases and groups loop over or use map function

phases = unique(set_1$points$TimePhase)

groups = unique(set_1$points$Role)

col_list_2 = c("red",#phase 1
              "blue", #phase 2
              "green") #phase 3

traces = c(2:10)

x = plot(set_1, title = "Group Means by Phase")
count = 1
for (i in 1:length(phases)){
  for (j in 1:length(groups)){
    points = set_1$points %>% filter(TimePhase == phases[i], Role == groups[j])
    x = x |>
      units(
        points = points,
        point_position_multiplier = point_position_multiplier,
        points_color = col_list_2[i],
        show_mean = TRUE, show_points = F, with_ci = FALSE
        ) |>
      nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = paste0(phases[i]," - ",groups[j]), traces = traces[count])
    count = count + 1
  }
}

x

NA

Statistical Tests:

In these data, the observations are nested within participants, so we need to check if this nesting has a significant effect. We do this by constructing a confidence interval around the intraclass correlation coefficient (https://en.wikipedia.org/wiki/Intraclass_correlation). If the interval contains zero, nesting is not significant.

reg_dat = set_1$points %>% filter(ENA_DIRECTION == "response")

ICC::ICCest(Speaker,SVD1,reg_dat) #not significant

Warning: 'x' has been coerced to a factor

$ICC
[1] 0.05242903

$LowerCI
[1] -0.2098135

$UpperCI
[1] 0.5007296

$N
[1] 10

$k
[1] 3.272727

$varw
[1] 0.3002277

$vara
[1] 0.01661158

ICC::ICCest(Speaker,SVD2,reg_dat) #not significant

Warning: 'x' has been coerced to a factor

$ICC
[1] -0.1818112

$LowerCI
[1] -0.3333495

$UpperCI
[1] 0.1964787

$N
[1] 10

$k
[1] 3.272727

$varw
[1] 0.07137526

$vara
[1] -0.01098045

Nesting is not significant so we proceed with an OLS regression. Below we set up two models for each dimension and test whether including interaction effects in the analysis significantly improves the modes.

mod_x = lm(SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat)
mod_y = lm(SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat)

mod_x_2 = lm(SVD1 ~ as.factor(TimePhase)*Role, data = reg_dat)
anova(mod_x,mod_x_2)

Analysis of Variance Table

Model 1: SVD1 ~ as.factor(TimePhase) + Role
Model 2: SVD1 ~ as.factor(TimePhase) * Role
  Res.Df    RSS Df Sum of Sq      F  Pr(>F)  
1     28 2.6341                              
2     24 1.7559  4   0.87815 3.0006 0.03851 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

mod_y_2 = lm(SVD2 ~ as.factor(TimePhase)*Role, data = reg_dat)
anova(mod_y,mod_y_2)

Analysis of Variance Table

Model 1: SVD2 ~ as.factor(TimePhase) + Role
Model 2: SVD2 ~ as.factor(TimePhase) * Role
  Res.Df    RSS Df Sum of Sq      F   Pr(>F)    
1     28 1.1117                                 
2     24 0.2156  4   0.89614 24.939 3.02e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Interactions are significant, so we should use the interaction models. Let’s view those models

summary(mod_x_2)

Call:
lm(formula = SVD1 ~ as.factor(TimePhase) * Role, data = reg_dat)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.85356 -0.01457  0.00790  0.08624  0.32816 

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)   
(Intercept)                                -0.1652     0.1352  -1.221  0.23385   
as.factor(TimePhase)2                       0.7013     0.2066   3.395  0.00239 **
as.factor(TimePhase)3                       0.6249     0.2066   3.025  0.00585 **
RoleProductDesigner                        -0.5044     0.1913  -2.637  0.01443 * 
RoleServiceDesigner                        -0.6041     0.1913  -3.158  0.00425 **
as.factor(TimePhase)2:RoleProductDesigner   0.1937     0.2815   0.688  0.49803   
as.factor(TimePhase)3:RoleProductDesigner   0.6236     0.2922   2.134  0.04324 * 
as.factor(TimePhase)2:RoleServiceDesigner  -0.2347     0.2815  -0.834  0.41272   
as.factor(TimePhase)3:RoleServiceDesigner   0.6444     0.2815   2.289  0.03117 * 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.2705 on 24 degrees of freedom
Multiple R-squared:  0.8261,    Adjusted R-squared:  0.7681 
F-statistic: 14.25 on 8 and 24 DF,  p-value: 2.062e-07

summary(mod_y_2)

Call:
lm(formula = SVD2 ~ as.factor(TimePhase) * Role, data = reg_dat)

Residuals:
      Min        1Q    Median        3Q       Max 
-0.196826 -0.035322 -0.001122  0.020930  0.273378 

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                               -0.17794    0.04739  -3.755 0.000977 ***
as.factor(TimePhase)2                     -0.08890    0.07239  -1.228 0.231326    
as.factor(TimePhase)3                      0.62627    0.07239   8.651 7.68e-09 ***
RoleProductDesigner                        0.24264    0.06702   3.620 0.001366 ** 
RoleServiceDesigner                        0.29823    0.06702   4.450 0.000168 ***
as.factor(TimePhase)2:RoleProductDesigner  0.06701    0.09865   0.679 0.503477    
as.factor(TimePhase)3:RoleProductDesigner -0.35748    0.10237  -3.492 0.001880 ** 
as.factor(TimePhase)2:RoleServiceDesigner -0.27608    0.09865  -2.799 0.009962 ** 
as.factor(TimePhase)3:RoleServiceDesigner -0.93799    0.09865  -9.508 1.30e-09 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.09478 on 24 degrees of freedom
Multiple R-squared:   0.89, Adjusted R-squared:  0.8534 
F-statistic: 24.28 on 8 and 24 DF,  p-value: 1.032e-09

Viewing interactions

emmip(mod_x_2, Role ~ TimePhase)

emmip(mod_y_2, Role ~ TimePhase)

When interactions are introduced, interpreting regression coefficients becomes more difficult, so we calculated the marginal means and various contrasts (i.e., comparisons between groups of interest) to address our research questions.

RQ1 (X dimension - Design Actions): Here we want to compare the “main effect” of Role on the x dimension—i.e., are there differences between the roles on the X dimension when averaging over TimePhases.

emm_x1 = emmeans(mod_x_2, specs = pairwise ~ Role, weights = "proportional")

NOTE: Results may be misleading due to involvement in interactions

print(emm_x1$contrasts)

 contrast                          estimate    SE df t.ratio p.value
 Engineer - ProductDesigner           0.251 0.119 24   2.116  0.1077
 Engineer - ServiceDesigner           0.487 0.116 24   4.192  0.0009
 ProductDesigner - ServiceDesigner    0.236 0.113 24   2.085  0.1142

Results are averaged over the levels of: TimePhase 
P value adjustment: tukey method for comparing a family of 3 estimates 

Only significant difference on the X is between Engineer and Service designer. Let’s calculate the effect size (difference in standard deviations) of that difference. We will use Cohen’s d.

cohensd = function(diff_,g1,g2){
  diff_/(sqrt((sd(g1)^2 + sd(g2)^2)/2))
}

contrasts = as_tibble(emm_x1$contrasts)
diff_ = contrasts$estimate[2]
g1 = reg_dat %>% filter(Role == "Engineer") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ServiceDesigner") %>% select(SVD1)

cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)

[1] 0.9779281

RQ1 (Y dimension - Design Actions): Here we want to compare the “main effect” of Role on the y dimension—i.e., are there differences between the roles on the X dimension when averaging over TimePhases.

emm_y1 = emmeans(mod_y_2, specs = pairwise ~ Role, weights = "proportional")

NOTE: Results may be misleading due to involvement in interactions

print(emm_y1$contrasts)

 contrast                          estimate     SE df t.ratio p.value
 Engineer - ProductDesigner          -0.157 0.0415 24  -3.771  0.0026
 Engineer - ServiceDesigner           0.078 0.0407 24   1.917  0.1556
 ProductDesigner - ServiceDesigner    0.235 0.0397 24   5.915  <.0001

Results are averaged over the levels of: TimePhase 
P value adjustment: tukey method for comparing a family of 3 estimates 

Significant differences between engineer/product designer and productdesigner/servicedesigner. Let’s get the effect sizes:

contrasts = as_tibble(emm_y1$contrasts)

diff_evp = contrasts$estimate[1]
g1 = reg_dat %>% filter(Role == "Engineer") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ProductDesigner") %>% select(SVD1)

d_evp = cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)

diff_pvs = contrasts$estimate[3]
g1 = reg_dat %>% filter(Role == "ProductDesigner") %>% select(SVD1)
g2 = reg_dat %>% filter(Role == "ServiceDesigner") %>% select(SVD1)

d_pvs = cohensd(diff_ = diff_,g1 = g1$SVD1,g2 = g2$SVD1)

d_evp

[1] 0.9137988

d_pvs

[1] 0.8072687

Marginal means - x RQ2

Here we use marginal means to test for differences between timephases within groups and groups within timephase and

library(emmeans)

emm_x1a = emmeans(mod_x_2, specs = pairwise ~ TimePhase|Role, weights = "proportional") #timephases w/n groups
emm_x1a = as_tibble(emm_x1a$contrasts)

emm_x2 = emmeans(mod_x_2, specs = pairwise ~ Role|TimePhase, weights = "proportional") #groups w/n timephases
emm_x2 = as_tibble(emm_x2$contrasts)

#see here: https://cran.r-project.org/web/packages/emmeans/vignettes/interactions.html

emm_x1a = emm_x1a %>% filter(p.value < 0.05)
emm_x2 = emm_x2 %>% filter(p.value < 0.05)

emm_x1a

emm_x2

NA

Several contrasts are significant. Let’s calculate the effect sizes

print(list(t1_evp,t1_evs,t2_evs,t2_pvs))

[[1]]
[1] 1.826995

[[2]]
[1] 2.233997

[[3]]
[1] 4.14715

[[4]]
[1] 1.169235

Some of these are too large to be trusted. This is happening because some individuals are highly similar in the data. We will have to acknowledge this as it is likely due to small sample size

Marginal means -Y RQ2: Let’s repeate what we did for X here

emm_y1 = emmeans(mod_y_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_y1 = as_tibble(emm_y1$contrasts)

emm_y2 = emmeans(mod_y_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_y2 = as_tibble(emm_y2$contrasts)

emm_y1 = emm_y1 %>% filter(p.value < 0.05)
emm_y2 = emm_y2 %>% filter(p.value < 0.05)

emm_y1

emm_y2

Several contrasts are significant. Let’s calculate the effect sizes

print(list(t1_ep_y,t1_es_y,t2_ep_y,t2_ps_y,t3_es_y,t3_ps_y))

[[1]]
[1] -2.082427

[[2]]
[1] -2.734181

[[3]]
[1] -13.53816

[[4]]
[1] 2.197361

[[5]]
[1] 8.972244

[[6]]
[1] 8.239644

Model diagnostics

library(performance)
library(see)

check_model(mod_x_2)

Model has interaction terms. VIFs might be inflated.
  You may check multicollinearity among predictors of a model without interaction terms.

check_model(mod_y_2)

Model has interaction terms. VIFs might be inflated.
  You may check multicollinearity among predictors of a model without interaction terms.

Things look ok for the most part

Plot network subtractions (replace TimePhase and Role for different combinations)

mean1 = set_1$points %>% dplyr::filter(TimePhase == 1, Role == "Engineer") #update here
mean2 = set_1$points %>% dplyr::filter(TimePhase == 1, Role == "ServiceDesigner") #update here

mean1_net = set_1$line.weights %>% dplyr::filter(TimePhase == 1, Role == "Engineer") %>% colMeans()
mean2_net = set_1$line.weights %>% dplyr::filter(TimePhase == 1, Role == "ServiceDesigner") %>% colMeans() #update here

plot(set_1, title = "Network subtraction") |>
  units(
    points = mean1, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = mean2, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = mean1_net - mean2_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))

Model 2: Epistemic Agency Codes

Plot 1: Mean Network

node_size_multiplier = 0.3 # scale up or down node sizes
node_position_multiplier = 1 # zoom in or out node positions
point_position_multiplier = 1 # zoom in or out the point positions
edge_arrow_saturation_multiplier = 1.5 # adjust the chevron color lighter or darker
edge_size_multiplier = 2 # scale up or down edge sizes

plot(set_2, title = "Overall Mean--Agency") |>
  units(
    points= set_2$points, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("black"),
    show_mean = TRUE, show_points = TRUE, with_ci = TRUE) |>
  edges(
    weights =set_2$line.weights,
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("black")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    node_labels = TRUE, # change this to FALSE can remove node labels in case you want to add them back in a nicer font or size for your presentations or publications
    self_connection_color = c("black"))

Plot 1b: Overall Group Comparisons

plot(set_2, title = "Groups") |>
  units(
    points= set_2$points$Role$ProductDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("red"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_2$points$Role$Engineer, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("blue"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
   units(
    points= set_2$points$Role$ServiceDesigner, 
    point_position_multiplier = point_position_multiplier,
    points_color = c("green"),
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = "Product Designer", traces = traces[1]) |>
      plotly::style(name = "Engineer", traces = traces[2]) |>
      plotly::style(name = "Service Designer", traces = traces[3])
  

Plot 1c: Group subtractions

prod_pts = set_2$points$Role$ProductDesigner 
eng_pts = set_2$points$Role$Engineer
serv_pts = set_2$points$Role$ServiceDesigner

prod_mean_net = set_2$line.weights %>% dplyr::filter(Role == "ProductDesigner") %>% colMeans()
eng_mean_net = set_2$line.weights %>% dplyr::filter(Role == "Engineer") %>% colMeans()
serv_mean_net = set_2$line.weights %>% dplyr::filter(Role == "ServiceDesigner") %>% colMeans()

plot(set_2, title = "Product Designer vs Engineer") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))

plot(set_2, title = "Product Designer vs ServiceDesigner") |>
  units(
    points = prod_pts, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points =serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = prod_mean_net - serv_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","green")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","green"))

plot(set_2, title = "Service Designer vs Engineer") |>
  units(
    points = serv_pts, 
    points_color = "green",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = eng_pts, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = serv_mean_net - eng_mean_net, # optional multiplier to adjust for readability
    edge_size_multiplier = 1,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("green","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("green","blue"))

Plot 2: All Means

#get list of time phases and groups loop over or use map function

phases = unique(set_2$points$TimePhase)

groups = unique(set_2$points$Role)


col_list_2 = c("red",
              "blue",
              "green")

traces = c(2:10)

x_agency = plot(set_2, title = "Group Means by Phase -- Agency")
count = 1

for (i in 1:length(phases)){
  for (j in 1:length(groups)){
    points = set_2$points %>% filter(TimePhase == phases[i], Role == groups[j])
    #print(head(points))
    x_agency = x_agency |>
      units(
        points = points,
        point_position_multiplier = point_position_multiplier,
        points_color = col_list_2[i],
        show_mean = TRUE, show_points = F, with_ci = FALSE
        ) |>
      nodes(node_size_multiplier = 0.3,
            node_position_multiplier = node_position_multiplier,
            self_connection_color = c("black"))|>
      plotly::layout(showlegend = TRUE, legend = list(x = 100, y = 0.9)) |>
      plotly::style(name = paste0(phases[i]," - ",groups[j]), traces = traces[count])
    count = count + 1
    
  }
}

x_agency

Check nesting

reg_dat_agency = set_2$points %>% filter(ENA_DIRECTION == "response")

ICC::ICCest(Speaker,SVD1,reg_dat_agency) #not significant

ICC::ICCest(Speaker,SVD2,reg_dat_agency) #significant we will need to include a variable for individual in the model

mod_x_agency = lm(SVD1 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)
mod_y_agency = lm(SVD2 ~ as.factor(TimePhase) + Role, data = reg_dat_agency)


summary(mod_x_agency)
summary(mod_y_agency)

check for interactions: (interactions are significant so we should use those modes)

mod_x_agency_2 = lm(SVD1 ~ as.factor(TimePhase) + Role + as.factor(TimePhase)*Role, data = reg_dat_agency)
mod_y_agency_2 = lm(SVD2 ~ as.factor(TimePhase) + Role + as.factor(TimePhase)*Role, data = reg_dat_agency)

summary(mod_x_agency_2)
summary(mod_y_agency_2)

anova(mod_x_agency,mod_x_agency_2)
anova(mod_y_agency,mod_y_agency_2)

#Need to run mixed model for Y ##model is singular so stick with OLS

library(lmerTest)

mod_y_agency_mixed = lmerTest::lmer(SVD2 ~ 1 + (1|Speaker) + as.factor(TimePhase)*Role ,data = reg_dat_agency)

summary(mod_y_agency_mixed)

Marginal means x - RQ1

emm_ag1 = emmeans(mod_x_agency_2, specs = pairwise ~ Role, weights = "proportional")
print(emm_ag1$emmeans)
print(emm_ag1$contrasts)

Viewing interactions

emmip(mod_x_agency_2, Role ~ TimePhase)

Marginal means - x RQ2

emm_ag1 = emmeans(mod_x_agency_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_ag1$contrasts

emm_ag2 = emmeans(mod_x_agency_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_ag2$contrasts

Marginal means y - RQ1

emm_ag_y1 = emmeans(mod_y_agency_2, specs = pairwise ~ Role, weights = "proportional")
print(emm_ag_y1$emmeans)
print(emm_ag_y1$contrasts)

Viewing interactions

emmip(mod_y_agency_2, Role ~ TimePhase)

Marginal means - y RQ2

emm_ag_y1 = emmeans(mod_y_agency_2, specs = pairwise ~ TimePhase|Role, weights = "proportional")
emm_ag_y1$contrasts

emm_ag_y2 = emmeans(mod_y_agency_2, specs = pairwise ~ Role|TimePhase, weights = "proportional")
emm_ag_y2$contrasts

Calculate eff sizes

#to do

Model diagnostics

check_model(mod_x_agency_2)
check_model(mod_y_agency_2)

Plot network subtractions

mean1 = set_2$points %>% dplyr::filter(TimePhase == 1, Role == "Engineer")
mean2 = set_2$points %>% dplyr::filter(TimePhase == 1, Role == "ProductDesigner")

mean1_net = set_2$line.weights %>% dplyr::filter(TimePhase == 1, Role == "Engineer") %>% colMeans()
mean2_net = set_2$line.weights %>% dplyr::filter(TimePhase == 1, Role == "ProductDesigner") %>% colMeans()

plot(set_2, title = "Network subtraction") |>
  units(
    points = mean1, 
    points_color = "red",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  units(
    points = mean2, 
    points_color = "blue",
    show_mean = TRUE, show_points = FALSE, with_ci = TRUE) |>
  edges(
    weights = mean1_net - mean2_net, # optional multiplier to adjust for readability
    edge_size_multiplier = edge_size_multiplier,
    edge_arrow_saturation_multiplier = edge_arrow_saturation_multiplier,
    node_position_multiplier = node_position_multiplier,
    edge_color = c("red","blue")) |>
  nodes(
    node_size_multiplier = node_size_multiplier,
    node_position_multiplier = node_position_multiplier,
    self_connection_color = c("red","blue"))